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The Fermi liquid- Wigner crystal transition in a two dimensional electronic system is revisited with 
a focus on the nature of the fixed node approximation done in quantum Monte Carlo calculations. 
Recently, we proposed (Phys. Rev. Lett. 94, 046801 (2005) ) that for intermediate densities, a 
hybrid phase (with the symmetry of the crystal but otherwise liquid like properties) is more stable 
than both the liquid and the crystal phase. Here we confirm this result both in the thermodynamic 
and continuum limit. The liquid-hybrid transition takes place at r* = 31.5 ±0.5. We find that 
the stability of the hybrid phase with respect to the crystal one is tightly linked to its delocalized 
nature. We discuss the implications of our results for various transition scenarii (quantum hexatic 
phase, supersolid, multiple exchange, microcmulsions) proposed in the literature. 

PACS numbers: 



I. INTRODUCTION 

The competition between electrostatic and kinetic en- 
ergy in a two dimensional electron gas is a problem which 
is simple to formulate (what is the phase diagram of N 
electrons on a surface S at zero temperature) yet diffi- 
cult to tackle. One reason for this difficulty lies in the 
smallness of the difference in energy between different 
phases (The Fermi liquid which is stable at high den- 
sity when the kinetic energy dominates^ and the elec- 
trostatically favored (Wigner) triangular crystal, stable 
at low densitj^). In the region of interest in this study, 
20 < r s < 80 (the dimensionless parameter r s oc \J S/N 
controls the ratio of the electrostatic energy over the ki- 
netic energy) where there is a strong competition be- 
tween the two kinds of energy, the difference of energy 
between the two phases is only of the order of 0.1% 
of the total energy so that very accurate methods were 
needed to study this problem. To illustrate this point, 
we plot in the Fig. 2] the energy of these two phases as a 
function of r s = m* e 2 / '(h 2 t^firn) (e being the electronic 
charge, e the dielectric constant, m* the effective mass 
and n = N/S the electronic density). At large distances, 
the physics is entirely controlled by electrostatic, and im- 
portant amount of energy are involved. Once we add a 
positive background to the system and make it globally 
neutral, we arrive at the energies plotted in the upper 
panel of Fig. ^ on the scale of the plot the liquid and 
crystal are completely indistinguishable. If we remove 
the Madelung energy (excess electrostatic energy of the 
uniform background with respect to the crystal) , we still 
cannot distinguish between the two phases (middle panel 
of Fig.^). It is only after we have removed the zero point 
fluctuation energy of the crystal that the curves become 
different and show a crossing at r s ss 37 where the Wigner 
crystallization has been believed to occur 3,4 (lower panel 
of Fig. Pi. 

In a seminal article in 1989, Tanatar and Ceperlej* 3 " 
used a fixed node quantum Monte Carlo (FN-QMC) tech- 
nique^ to locate the critical value r s ss 37 ± 5 where the 




FIG. 1: Quantum Monte-Carlo energy as a function of r s for 
a system of 56 electrons in 210 x 208 sites for the liquid phase 
(squares), solid phase (diamonds) and hybrid phase (full cir- 
cles). Upper panel: total energy, the three phases are indis- 
tinguishable. Middle panel: total energy minus the Madelung 
(electrostatic) energy — c\/r a with c\ — 2.2122. Lower panel: 
total energy minus the Madelung energy and the energy of 
the crystal's phonons c 1 / 2 ^s~ 3/ ' 2 with C1/2 = 1.624. The en- 
ergies are given in mRy. For holes in a GaAs heterostruc- 
ture, a Rydberg corresponds roughly to 350 Kelvin (e ~ 13eo , 
m* ~ 0.38m e ). 



quantum melting of the Wigner crystal occurs. Their 
work was followed by more precise numerics^ and a bet- 
ter description of the liquid phase 7,8 that included back- 
flow corrections. On the other hand there were indica- 
tions that this scenario of a (first order) direct transition 
between the Wigner crystal and the Fermi liquid phase 
might miss some of the physics. For instance, the quan- 
tum melting of the bosonic Wigner crysta&iS is found 
(by similar QMC calculations) to occur at r s w 60 so that 
for 37 < r s < 60 the fermionic statistics (from which all 
the difficulties of the QMC calculations come, as we shall 
see) is crucial to stabilize the crystal. Also, the classical 
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meltingiijiSjiSjU ag a f unc ti n of temperature does not 
occur in a one step process. The system first looses its 
translational order but retains some orientational order 
(hexatic phased) while at a higher temperature, any or- 
der disappears. The possibility of an intermediate quan- 
tum hexatic phase was put forward in^A^. In fact, it 
was recently argued on rather general grounds that a 
direct first order transition is simply impossible in this 
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but should occur through a serie of interme- 



diate phases (including bubbles or stripes of one phase in 
the other). The possibility of a supersolid phase analo- 
gous to the one proposed by Andreev and LifshitaiS has 
also been considered22i2i*2£. All these indications— lead 
us recently to revisit the melting of the Wigner crystal 
in the framework of the FN-QMC calculations. We in- 
troduced a third phase called hybrid phase (as it had the 
symmetry of the crystal but is otherwise delocalized) and 
showed that it was more stable than both the liquid and 
crystal phase (Full circles in the lower panel of Fig. ^| in 
the intermediate density ranged giving a strong support 
to the more exotic scenarios. 

To understand the status of our result however, it is of 
prime importance to understand the nature of the (fixed 
node, FN) approximation involved in the FN-QMC tech- 
nique. In practice, the FN-QMC algorithm is fed with 
a wave-function called the guiding wave function (GWF) 
that has to be given explicitly, and that should be as 
close as possible to the ground state of the system. The 
FN-QMC algorithm projects the GWF onto the (true) 
ground state of the system. However, to avoid the no- 
torious sign problem that arises in QMC simulations in- 
volving fcrmions, one approximation is introduced: the 
projection is done with the constraint that the sign of 
the wave function (nodal structure) remains unchanged 
at every point of the Hilbert space. The method gives 
the best wave function for a given structure of the nodes 
of the GWF and is in this sense variational— The phys- 
ical meaning of the FN approximation is not completely 
obvious. Implicit in the interpretation given above of the 
FN-QMC results is that a GWF is associated to a phase 
(crystal, liquid or hybrid) the stability of which can there- 
fore be studied. However, this paradigm should not be 
taken too literally and will be revisited in this article. 

After presenting our model and the FN-QMC method 
in section II, we perform FN-QMC calculations on a small 
system of 4 electrons (section III) . As this system can be 
studied exactly, it is an interesting tool to study the na- 
ture of the FN approximation. We find that in addition 
to the FN results at large imaginary (projection) time, 
important information is embedded in the evolution of 
the results from the variational calculation to the full FN- 
QMC results. In section IV, we explain the construction 
of the hybrid GWF, and show that its nodal structure 
corresponds to delocalized waves. Section V is devoted 
to a characterization of the physics associated with the 
hybrid GWF. This is done through a systematic study 
of various physical quantities, and in particular of their 
evolution between the variational and FN-QMC calcula- 



tions. We find that the success of the hybrid GWF with 
respect to the crystal one is closely linked to its delocal- 
ized nature. Section VI contains a detailed discussion of 
various technical aspects (thermodynamic limit, lattice 
effects, mixed and unbiased estimates,...) The chief re- 
sult of section VI is a precise determination of the critical 
value r*=31.5±0.5. at which the liquid- hybrid transi- 
tion takes place (in both the continuum and thermody- 
namic limit). In the discussion section VII, we discuss 
the implications of our results to the scenarii proposed in 
the literature. 



II. MODEL AND METHOD. 

A. System under consideration 

The system under consideration consists of N spin- 
less electrons on a rectangular L x x L y grid with nearest 
neighbor hopping and long range Coulomb repulsion. To 
avoid strong finite size effects on the electrostatic energy, 
the system is repeated periodically and fills the whole two 
dimensional plane. In practice, we use periodic bound- 
ary conditions for the hopping terms and the effective 
two-body interaction is obtained from the bare Coulomb 
interaction using the Ewald summation technique. The 
system Hamiltonian reads** 

H = -tJ2 4 c r' + r'hrnp + A. (1) 

(r,r') r^r' 

where the operator ct (c?) creates (destroys) an electron 
on point r with the standard anticommutation relation 
rules, the sum J^(r f) ' IS done on the nearest neighbor 
points on the grid and t is the corresponding hopping 
amplitude. The density operator reads n? = etc?. U is 
the effective strength of the two body interaction V(r) 
which reads, 
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In the previous equation, k c is a (irrelevant) cut off. The 
vector L takes discrete values L = {n x L xl n y L y ) with 
n. x and n y integer numbers. The vector K also takes 
discrete values, K = (j^-n x , ff-%) and (n x ,n y ) ^ (0,0). 
The complementary error function is defined as Erfc(r) = 
fx° e_t dt. In order to assure electrostatic neutrality 
we add a positive continuum background 



X/N = At + UV(6) - 2Uuy/n/k c - 2Uk c 



(3) 
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where v = L L is the average number of electrons per 
site and V(r) = V(r) with the restriction that the sum 
over L does not includes the null vector. 

The presence of the grid can be understood either 
as a discretization of the continuum problem (and the 
nearest neighbor hopping corresponds to the discretized 
Laplacian) or as a tight-binding approach to two di- 
mensional electron systems where each site corresponds 
to an atomic orbital. In conventional 2d systems (say 
GaAs/GaAlAs heterostructures) , densities range from 
n s w 10 12 cm -2 down to n s « 10 10 cm~ 2 . Since the 
distance between nearest atoms is of the order of a few 
Angstroms, this would lead to v ~ 10~ 4 to v ~ 10~ 2 elec- 
trons per site in our tight-biding picture. In this study, 
we will study systems with v = 1/56, v = 1/224 and 
v = 1/780 close to realistic values. 

The presence of the jeHium (which is merely a con- 
stant term and thus cannot affect the physics) allows 
us to make a quantitative contact with the literature 
in the continuum limit as v — ► 0. In this limit, the 
physics depends only on the r s parameter which reads, 
r s = U I (2t^fiu?) while the Rydberg unit of energy is 
Ry = U 2 1 (At). In the following, unless specifically stated, 
we shall measure all energies E in unit of 2-nNvt (i.e. 
energy per particle in unit of the Fermi energy of the 
non-interacting problem.) With this normalization, the 
energy of the system at r s = in the thermodynamic 
N 3> 1 and continuum v <C 1 limit is E = 1 so that 
results in Rydberg can be obtained by multiplying the 
energies (in unit of 2tt Nut) by 2/r 2 . 

We have added a very small disorder ^2pVpnp in or- 
der to lift the degeneracies of the non interacting prob- 
lem, vp are independent and uniformly distributed inside 
[—W/2,+W/2]. We choose W — 10~ 3 corresponding to 
an extremely large ratio l/Xp = 96ut 2 /W 2 ~ 10 6 of the 
mean free path I over Fermi wave length Xp. We explic- 
itly checked that our result are insensitive to the presence 
of this disorder. 

The grid has been chosen to accommodate a (almost) 
triangular Wigner crystal without distortion. Hence we 
use systems of N = Q x R electrons (R (R even) lines of 
Q particles) in Qd x x Rd y sites. The triangular crystal 
requires d y /d x = V3/2 « 0.866. In this work, we focus 
on d y /d x = 14/16 w 0.875 and d y /d x = 26/30 » 0.866 
for which we find that the distortion is negligible. 



B. The Green Function Monte-Carlo method in 
the Fixed Node approximation 

We aim to sample the ground state |$o) 01 H and its 
corresponding ground state energy Eq with the Green 
Function Monte-Carlo technique. The idea behind this 
method is to project an initial (variational) guiding wave- 
function \^>g) on the exact ground state |^o) by applying 
the operator e~ H/3 in a stochastic way. In practice, one 
works in the many particle basis R — (r*i, r2---rjv) and 



applies the function 

G W r = Sr-r - tVg(R') [Hwr - u8 RIR ] y G \R) (4) 

where *S?g(R) — (R\^g), t is a (small) time step and cu 
is an unimportant offset of the energies (that should be 
roughly set to the ground state energy of the system, here 
we took u)/(2irvt) = l — r a ). In the absence of the ^g(R), 
G is a discretized version of the operator e~ H P for a small 
(imaginary) time step r. Upon applying n = /3/t times 
G on a vector ^ 2 G (R), the ^g(R') on the left of Eq.Q are 
canceled by the (R) on the right so that one obtains 

I n = J2(G n )R'R^ 2 G (R) = <* G | e - (ff - w) " r |*G> (5) 
R'R 

from which one can extract the ground state energy (as 
I n+ l = e- [En -^ T I n when (3 = nr -> oo). 

The stochastic implementation of this scheme is based 
on the Green Function Monte Carlo for lattice Hamilto- 
nians introduced in Ref. |2(| to which we refer for more 
details. The algorithm to update the Slater determinants 
used in the calculation of ^/g(R)/^g(R') can be found 
m Ref. HI By sampling directly the time spent by the 
walkers at one point of the Hilbert space using the al- 
gorithm described in Ref 26| we can use arbitrary small 
time steps r without any loss in computing time and 
hence effectively work in continuum (imaginary) time. 
Instead of using the standard branching technique, the 
control of the walkers population is done using a fixed 
number of walkers and the reconfiguration algorithm in- 
troduced by Sorella 28 . This algorithm allows to avoid the 
bias introduced in the branching technique by artificially 
controlling the walker population. Quantum averages of 
physical quantities (...) are calculated using the forward 
walking technique 2 ^, and hence do not suffer from the 
bias of mixed estimates. 

So far the scheme is essentially exact. However, it 
suffers from the usual "sign problem", the sign of Grir 
fluctuates so that the statistical accuracy decreases ex- 
ponentially with /3, and it is of little practical use. One 
way out of the sign problem is the Fixed-Node approx- 
imation where one forbids the sign of the wave func- 
tion to change upon applying G (hence the name: the 
nodal surface where the wave-function changes of sign 
remains the same than the one of ^g(R))- The practi- 
cal implementation of the fixed-node approximation on 
a grid is done 25 by replacing H by an effective Hamil- 
tonian _ff FN that depends on the GWF. H™ R is equal 
to H R ' R when G R ' R > 0. When G R ' R < 0, the link 
is cut H R r R = and is replaced by an effective poten- 
tial H™ = Hrr + J2 R < 6(-G R .R)* G (ff)H R , R *G t (R) 
where 0(x) is the Heaviside function (this corresponds 
to replace the value of the wave-function on those sites 
R! by the surmise ^g(R'))- The fixed node approxima- 
tion can be thought as a "supervariational" technique 
where the amplitude of the wave-function is optimized 
at every point of the Hilbert space while its sign remain 
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FIG. 2: (color online). Typical curve of £fn(/3) as a function 
of p. System of 32 particles in 64 x 112 sites at r s = 40 for 
a liquid GWF with A = 4.5. Here, 8.10 walkers were used, 
uj — — 39 and r = 10~ /t. Inset: same curve as a function of 
1//3. The dashed lines are linear fits to guide the eye. 



fixed. It can be proved indeed that the energies -Efn cal- 
culated with i?FN are larger than the true ground state 
energy Eq but smaller than the variational energy asso- 
ciated with the guiding wave-function2^. The fixed node 
approximation in a lattice looks a priori more drastic 
than in the continuum since we do impose the ratio of 
the wave-function across the nodal surface. However, for 
f <gc 1 the fraction of "nodal" sites goes to zero and the 
technique becomes equivalent to the fixed node diffusive 
Monte-Carlo used in the continuum 3 . 

As an illustration, we present in Fig. a typical trace 
obtained for 32 particles at r s — 40. Efn(/? = 0) corre- 
sponds to the variational energy c\H\^> a) / g\^ g) ■ 
After an initial rapid decrease, the FN energy (per parti- 
cle) Efn(P) decreases slowly as 1/(3 (see inset) and then 
saturates above an imaginary time (3 sa t cx l/(yt). The 
energy is estimated by further averaging the result over 
(3 for j3 > 0sat- The typical achieved accuracy allows 
us to determine £fn(/3) with a precision better than 
±0.001. This remarquable precision (here the relative 
accuracy is ~ 10 -5 ) should be contrasted with the fact 
that important changes in the physics can lead to very 
tiny changes in energy. For instance the condensation 
energy in a superconductor is only a very small fraction 
~ (A/Ep) 2 ~ 10 -6 (A superconducting gap, Ep Fermi 
energy) of the total energy (only electrons near the Fermi 
surface form Cooper pairs) . The formation of a Wigner 
crystal (localized in real space) from a Fermi liquid (lo- 
calized in momentum space) corresponds however to a 
complete reorganization of the system and the changes 
in energy, though small ~ 0.01, can be measured with 
the GFMC technique. 



C. Guiding wave functions. 

Central in the Fixed-Node technique is the choice of 
the GWFs used in the calculations. Here, the guiding 
wave-functions have the general form of a Slater deter- 
minant multiplied by a Jastrow function, 



9 G (R) =Det [fcirj)] x JJjfli 



(6) 



The Jastrow part takes Coulomb interaction into account 
by introducing correlations between electrons. Since it 
has no nodes, the FN results should not depend on its 
particular form (as we explicitly checked). We use mod- 
ified Yukawa functions^, 



J(r) = exp 



aA(r s ) 



(1 



-B(r s )r/ a) 



(7) 



where a = 1/^/ttv is the average distance between elec- 
trons. A(r s ) and B(r s ) are (optimized) variational pa- 
rameters. To avoid the Coulomb singularity when two 
electrons get close to each other , we impose the cusp con- 
dition that reads B = \/r s /A for the modified Yukawa. 

The Slater determinant of one-body wave functions, 
Det [<pi (fj)] enforces the antisymmetric nature of the 
fermionic wave function and is responsible for the nodal 
structure of the GWF. The GWF used in the literature 
for the study of the melting of the Wigner crystal 3 ^ are 
of two kind, adapted to the two limits of very low (large) 

• For the liquid state the GWF is constructed 

out of plane waves 4>iifj) cx e lki ' Tj with a well de- 
fined Fermi surface. \&i iq is the exact ground state 



at r s = 0. 

For the crystal GWF 

(piifj) (X e -(rj-Si) 2 /(2d 2 ) 



$ rrv , localized orbitals 
Here the 



cry j 

are used 



with i G {1 . . . N} stand for the positions of the 
electrons in the classical crystal. ^ C i y provides the 
exact ground state in the continuum at large r s 
with the variational parameter do cx a/rl^ 4 . This 
GWF captures the two leading terms of the large r s 
expansion of the energy (Madelung and zero point 
fluctuation energy) that reads, 

E = - Cl r s /2 + c 1/2 ^T s /2 + --- (8) 

with d = 2.2122 and c 1/2 = 1.624. 

In this article, we will use the two previously mentioned 
GWF, and introduce a third (hybrid) one that somehow 
interpolates between the crystal (real space) and the liq- 
uid (momentum space). 



III. THE SIGN PROBLEM AND THE FIXED 
NODE APPROXIMATION 

In this section, we come back to the sign problem, and 
investigate the nature of the approximation involved in 
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FIG. 3: Schematic view of Eq.(|5J. The circles symbolize sites 
7? in Hilbert space and the full (dashed) lines positive (nega- 
tive) off-diagonal matrix elements H RR i. The corresponding 
lattice is frustrated for fermions. 



the Fixed Node approximation. One way to understand 
the sign problem is to consider it as a frustration prob- 
lem in Hilbert space: we seek to solve the Schrodinger 
equation which formally reads in R space, 
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H RR ^ (Bf) = E ^ (R). 
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FIG. 4: Comparison between the Fixed Node £pn (upper 
panels) and variational -Evar (lower panels) energies as a 
function of do (left panels) and A (right panel) for a system of 
30 particles in 40 x 42 sites at r s = 35 using the * cry GWF. 
Left panels: A — 4.5. Right panels: do = 0.76. 



Two-body interactions and external potentials are diag- 
onal in the R space and Hrr can always be considered 
to be positive (by shifting H by a constant). The off- 
diagonal elements Hrr> come from kinetic energy. They 
are positive for bosons but are alternatively positive and 
negative for fermions as shown schematically in Fig. [3] 
The bosonic case does not suffer from any sign prob- 
lem and ^E , Q OSOn (_R) can be found efficiently. We now 
consider a fictitious Ising model where an fsing vari- 
able s(R) = ±1 is placed on each site 7? of the Hilbert 
space and is coupled to its neighbors by a coupling Hrri 
(which can be ferro or antiferromagnetic) . If this fic- 
titious model had no frustration, then its ground state 
sq(R) — ±1 could be trivially found and one can eas- 
ily verify that * crmion (i?) = s (R)^ oson (R) would be 
the exact ground state of the fermionic problem. Hence 
the sign problem arises from the frustrated nature of 
Eq.@. From this point of view, the Fixed-Node approx- 
imation consists in cutting some links of this fictitious 
fsing model so that it is no longer frustrated and ad- 
mits sq(R) = sgn(^/a(R)) as its ground state. The FN 
approximation does not have a simple physical meaning 
and thus seems difficult to control. 

In Fig. 01 we compare the FN results with the varia- 
tional (i.e. (3 = 0) results for a system of 30 particles at 
r s = 35. The calculations have been done for the ^cry 
GWF so that we can compare the influence of the two 
parameters do and A. While do enters in the definition 
of the Slater determinant and hence of the nodal struc- 
ture, A defines the Jastrow part which does not change 
sign and (in the limit v -C 1) the FN results should thus 
be insensitive to the choice of A. On the right panel of 
Fig. 0] we see that it is indeed the case. The function 
E-p-^(A) is flat at the precision of the calculation (±0.002 
here) while the variational energy Eyar(A) shows a pro- 
nounced minimum at A ss 4. On the left panels of Fig. 21 
we plot -Efn and -Evar as a function of do- One would 



expect the two curves to present a minimum at the same 
value of do (which is what is usually assumed, the FN cal- 
culations being done in most cases with parameters op- 
timized at the variational level). However, one observes 
that for do > 0.8 (variational minimum) the FN energy is 
flat upon increasing d - Hence, the nature of the approx- 
imation done in a fixed node calculation is not simply 
related to the original variational wave-function used for 
fixing the nodal surface, and the method can potentially 
capture more physics than originally thought. 

To get further insight, we now discuss the instructive 
example of 4 electrons in a 6 x 8 grid. This system is small 
enough (the size of the Hilbert space is ~ 2.10 5 ) so that 
its ground state can be found exactly using the Lancsos 
algorithm and compared with the result of a Fixed Node 
calculation. We introduce the density p(r) and density- 
density correlation function g(f) 



P(r) = - 
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(10) 
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which measures respectively the average number of elec- 
trons on site r and the probability to find an electron on 
site r knowing that there is one particle on site (both 
are normalized to 1). From these two quantities, we de- 
fine the participation ratio 



P2 = l/ 



L x L y J2p 2 (r) 



(12) 



which roughly measures the number of occupied sites in 
the ground states (p2 interpolates from 1 to v as r s goes 



from to oo ) and the "rigidity" parameter gi 



fj2 



L X Ly 
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The comparison between the exact result and the FN 
calculation is presented in Fig. [S] In the inset of Fig. 
we show the exact g(r) at r s — 20. g(f) has three 
well defined maxima at the positions of the classical 
Wigner crystal. Upon increasing r s , the system gets more 
and more crystal like and these maxima gets more pro- 
nounced which leads to an increase of 32 (lower panel) 
and a decrease of P2 (middle panel) . Figure |S] shows 
two regimes. For r s < 30, the error on the energy (up- 
per panel) increases linearly while the two other physical 
quantities P2 and gi are in very good agreement with the 
exact result. Hence the FN result, though calculated with 
a liquid GWF made of plane waves, is well able to repro- 
duce the ground state of a system which (at r s = 30) is 
really crystal like. Above r s > 30 the error saturates to 
0.007, P2 is still in good agreement with the exact result, 
but the error on gi starts to increase significantly. At 
r s = 70, the error on the energy is only 0.01% of the 
total energy (Eq rs 1 — r s ) while the error on gi is much 
larger ~ 33%. 

Hence the FN approximation provides a very good de- 
scription of the ground state up to a rather large value 
of r s (here around 30) even though the system is very 
crystal like and we started from a liquid GWF. It also 
provides a good quantitative value of the energy even for 
higher values of r s where the FN approximation no longer 
describes accurately the physics. Maybe even more im- 
portantly, Fig. gives us a clue on how to control the 
FN approximation: we find that when the FN approx- 
imation starts to fail (r a ~ 30), the difference between 
the FN results and the variational results starts to in- 
crease drastically. In addition, the FN results (here for 
pi and g-i) lies in between the variational and the ex- 
act results so that the fixed nodal surface can be viewed 
as a "wall" preventing the system to relax to the true 
ground state. Hence some information on the nature of 
the ground state can be extracted from the evolution of 
the physical quantities between the variational and the 
FN estimates (and not only from the study of the FN 
quantities alone) This idea will be put to application in 
section by looking how the physics put in the various 
GWF is amplified or washed out by the FN projection. 



IV. CONSTRUCTION OF AN HYBRID PHASE 

The melting point of the Wigner crystal can be viewed 
as the point above which the problem is well described 
in real space (crystal with particles at given positions) 
and below which the momentum space is to be used (liq- 
uid with well defined Fermi surface). Hence, the melting 
point itself is somehow the point of "maximum uncer- 
tainty" where the problem is equally badly described in 
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FIG. 5: (color online). Comparaison between the Fixed Node 
calculation with a liquid GWF and the exact result for a sys- 
tem of 4 particles in 8 x 6 sites. Upper panel: difference 
between the FN and exact energy. The lines are guide to the 
eye. Middle panel: participation ratio per site p2 for the exact 
(empty circles), FN (full diamonds) and variational (crosses) 
calculation. Lower panel: idem for the rigidity parameter 
(/2- Inset: Exact density-density correlation function g(r) at 
r s — 20. g(f) measures the probability to find a particle on 
site r knowing that one particle is on site indicated by the 
cross in the corner of the sample. The grayscale range from 
(black) to 3 (white). 



both momentum and real space representation. In the 
standard scenario of a (first order) direct transition be- 
tween the liquid and the crystal, this point is nothing else 
than the point where the energy of one state crosses the 
other one. Here, we consider another scenario which in- 
terpolates between momentum and real space. In order 
to do so, we introduced in Ref. a new GWF, called 
hybrid GWF aimed to provide such an interpolation. As 
we have seen, we find that this hybrid GWF provides the 
lowest energy in the intermediate region of r s . 

The construction of the hybrid GWF ^hyb is done in 
such a way that the resulting (f>i(fj) are the Bloch states 
of electrons in a triangular crystal. First, an effective 
one-body Hamiltonian H c g is constructed for an electron 
in a attractive periodic potential that has the symmetry 
of the classical (triangular) Wigner crystal, 



H c s — —t 



(r,r') 



r )np 



where the one-body potential is W(r) = YU=i Ui). 
The singularity of W(r) at the position of the classical 
crystal r — Uj has been removed by setting W(ui) = 
W(u.i + (1,0)) (we checked that other choices of W(r) 
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FIG. 6: Fixed Node energy (left panels) and variational en- 
ergy (right panels) at r a = 15 (upper panels) and r 3 — 35 
(lower panels) as a function of U* . The system contains 30 
particles in 40 x 42 sites and the calculations are done with 
the hybrid GWF. 



rower as U* increases. What is important here is that the 
point where this band of lowest energy separates from the 
rest of the level (see Fig[7| upper paner which is a zoom 
of the inset) corresponds exactly to the point where the 
energy has a drop (see Fig[7| lower paner). Above this 
point, Bloch theory tells us that the <fii(r) can be written 
in term of the Bloch waves of one band only: 



(15) 



where (r) is a function with the periodicity of the crys- 
tal and the momentum have their value within the first 
Brillouin zone of the triangular lattice, i.e. an hexagonal 
Fermi surface. In that sense, the hybrid GWF is liquid 
like (made of plane waves) with the triangular symmetry. 
However, since the only thing that enters in the GWF is 
actually the determinant, det[^i(r)], it is equivalent to 
use (instead of Eq. i|15|) ) linear combinations of the Bloch 
waves, leading to the so called Wannier functions: 



4>i(f) = F{r — Uj) with F(r) 



V(r>^- (16) 



give consistent results). In a second step, H e g is numer- 
ically diagonalized using the Lanczos algorithm. The N 
orbitals of lowest energy <pi(r) (1 < i < N) are then 
used to construct the Slater determinant. A priori U* is 
a variational parameter that allows for an interpolation 
between "Fi; q (at U* = 0) and localized orbitals leading 
to a GWF similar to * cry (at U* > 1). 

In FigEl we plot the energy as a function of U* for 
two values of r s below (r s — 15 upper panels) and above 
(r s = 35 lower panels) the expected liquid-crystal tran- 
sition (which is found around r s w 20 — 25 for this sys- 
tem of 30 particles in 40 x 42 sites). At r s — 15 the 
energy increases with U* indicating that the free plane 
wave solution provides the best nodal structure. The in- 
crease of the Fixed-Node energy (left panels) is however 
10 times smaller than the increase of the variational en- 
ergy and is only slightly above the statistical resolution. 
At r s = 35, the situation is completely changed, and the 
FN energy shows a sharp drop around a particular criti- 
cal value U* = U c « 0.2. This point corresponds, as we 
shall see, to the value of U* above which the GWF has 
the symmetry of the crystal. We find that the fixed node 
energy is completely flat above and below this threshold, 
indicating that the drop in energy is really associated to 
the change of symmetry. 

More insight on the significance of the critical value U c 
can be found by looking at the eigenenergies ej associated 
with the orbitals <j>i{r). The lowest values of these ener- 
gies 6i, which have no physical meaning by themselves, 
are plotted in the inset of Fig. [7| as a function of U* . 
At U* = they form one unique band with a parabolic 
dispersion a cx kf . At U* 3> 1, the periodic potential is 
at the origin of a more complex band structure, and the 
band of lowest energy (which contains exactly 7Y levels) 
detaches from the rest of the spectrum and becomes nar- 



Writing the <j>i(r) in the form Ea. (jrp|) . the hybrid GWF 
looks similar to the crystal one (for which F{r) takes the 
form of a gaussian function). The Wannier functions are 
not uniquely defined (since (r) is defined up to an ar- 
bitrary phase) and an important effort has been done in 
the literature (in the context of ab initio electronic struc- 
ture calculation mostly) to define the maximally localized 
Wannier functions. The fondamental result for us is22i2i 
that at U* = Uq, F{f) cannot decrease faster than al- 
gebrically with r and is hence a delocalized function over 
the sample. This is in sharp contrast with the crystal 
GWF for which the (f>i(r) are strongly localized around 
the classical positions Ui of the crystal. For U* Uq, the 
lowest band becomes narrower, and exponentially local- 
ized Wannier function F(f) ~ e~ hr can be constructed. 
However, the point of interest for us is U* w U c where 
the drop in energy occurs. In the rest of this paper, we 
will stick to this value where the Wannier functions are 
delocalized and hence completely different from the crys- 
tal GWF at a qualitative level. 



V. CHARACTERIZATION OF THE HYBRID 
PHASE 

As we have seen in the introduction (lower panel of 
Fig. nj, the hybrid GWF provides the lowest energy for 
a rather large domain of r s . Having established that, we 
now try to understand the physical origin of this success 
and the nature of the corresponding phase. From the con- 
struction of the hybrid GWF (Slater determinant made 
of a delocalized Wannier function) , we expect (at the vari- 
ational level) something someway between the crystal (it 
has its triangular symmetry) and the liquid (it is made 
of delocalized one-body functions) . While the crystal as- 
pect is straightforward, the characterization of the liquid 
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FIG. 7: (color online). Upper panel: lowest energy levels 
£i of the fictitious one-body problem H e a as a function of 
the variational parameter U* for a system of 72 particles in 
96 x 168 sites. The dashed lines correspond to the 72 th and 
73 th values of £j. They indicate the critical value U* ~ 0.12 
where the band of lowest energy splits from the rest of the 
fictitious spectrum. Inset: same thing for a broader range of 
U* . Lower panel: variational energy as a function of U* for 
2P 2 particles in 16P x 28P sites at r a = 40 with A = 4.5. 
P = 3 (empty diamonds) , P = 4 (empty squares) , P = 5 (full 
circles) and P = 6 (full diamonds) . 



aspect is a difficult task within the FN-QMC. Indeed, 
the FN-QMC method mainly provides access to quanti- 
ties (like p(f) or g(r)) which are diagonal in the real space 
representation, and hence not very likely to characterize 
the liquid behaviour. 

In Ref. |3 we studied some physical quantities for the 
three GWF with the implicit idea that these quantities 
would characterize the phase associated with each GWF. 
Here we follow a different approach. Following the in- 
sight we got in section III on a small system, we study 
how the physical quantities evolve between the varia- 
tional value Ovar = (*g|0|*g) and the FN results 
Ofn = (*o N l ( ^ ) l v ''o N )' as tn i s g lves the general trend 
toward the true ground state of the system. The differ- 
ent GWFs capture different physics, which can be either 
amplified or washed out by the FN projection procedure, 
giving hints on the true nature of the ground state. 

In the right panel of Fig|S|we plot the kinetic and inter- 
action energies as a function of r s . The variational values 
of the hybrid GWF are in the middle of those obtained 
with the liquid and crystal GWF. Upon applying the 
FN projector e~@ H on the hybrid (and crystal) GWF, 
the value of the kinetic energy decreases and converges 
to a value rather close to the variational kinetic energy 
of the liquid. It is important to note that these differ- 
ences in (kinetic or interaction) energy are about 5 to 10 
times bigger than the differences observed on the total 
energy. This decrease of the kinetic energy to an almost 
liquid value is a strong sign of the delocalized nature of 
the ground state. As r s increased however, this tendency 
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FIG. 8: Evolution from variational to FN results for various 
physical quantities. p m i n (r s ) = mhvp(r) (Upper Left panel), 
Pnmx(r s ) = max ? p(f) (Lower Left panel), E c (r a ) (Upper 
Right panel) and Ei n t(r 3 ) (Lower Right panel, the Madelung 
energy — c\r s j2 has been substracted) for a system of 56 elec- 
trons in 56 x 56 sites. Average done on 5 x 10 walkers. The 
empty symbols show the variational results while the full sym- 
bols show the FN walking estimates at — 16 for the hybrid 
(diamonds), crystal (square) and liquid (circles) GWF. The 
arrows indicate the evolution from the variational to the FN 
estimate. p m i n and p max are measured in percentage of the 
average density. 



also diminishes. This delocalization can also be seen indi- 
rectly on the electronic density p{r). An example of p{f) 
is shown in the upper panel of FigBJ At first sight the 
density looks crystal like (strong peaks forming a triangu- 
lar lattice) . However the height of the peaks (p ma , x lower 
left panel of Fig 0) and the depth of the valleys ( p m i n up- 
per left panel of Fig|8J) show not only that these peaks are 
not very pronounced (at r s = 32, the density at the peaks 
is only twice the average density while the "background" 
contains p m [ n w 35% of the electrons) but also that the 
contrast tends to diminish from the variational to the 
FN estimate (especially for the crystal GWF). Similar 
conclusions can be drawn from the density-density corre- 
lation function g(r) (A 3D plot is shown in Fig|SJ FigllOl 
gives a cross section). For r s — 35 (upper panel of Fig llTJ|) 
the peaks of both the hybrid and crystal results tend to 
be washed out after the FN projection At r s = 60 how- 
ever (lower panel of Fig llO|) the tendency is inversed and 
the hybrid result tends to get more localized when going 




FIG. 9: (Color online.) Density p(r) (Upper panel) and 
density-density correlations g(f) (Lower panel) at r s = 35 
for 56 electrons in 56 x 56 for the hybrid GWF. The density 
is plotted in percentage of the average density. Its minimum 
(maximum) correspond to 35% (230%) of the average density. 



from variational to the FN results (the same thing can 
be seen on p max , see lower left panel of FigjSJ. 

All the physical quantities discussed above indicate the 
same tendency: at the variational level the hybrid GWF 
lies somewhere between the liquid and the crystal results, 
and the FN results get even closer to the liquid for not 
too large r s (roughly r s < 60). At a basic level, the raw 
reason why the hybrid GWF allows to gain energy with 
respect to (the previously used) crystal GWF is therefore 
quite simple: electrons in the crystal GWF are far too 
much localized, and the hybrid GWF helps to restore a 
better balance between kinetic and electrostatic energy. 



VI. TECHNICAL ASPECTS 

In this section, we address various technical points. We 
start with giving some reference values of the energy that 
can be compared directly with the literature. We proceed 
with studying the thermodynamic (N — > oo) and contin- 
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FIG. 10: Cross section of g(r) with r — {x, 0) as a function 
of x for 56 electrons in 56 x 56. The empty (full) symbols 
are the variational (FN) results, respectively for the hybrid 
(circles), crystal (diamonds) and liquid (squares) GWF. The 
upper (lower) panel corresponds to r s = 35 (r s — 60). The 
arrows indicate the direction of the evolution from the varia- 
tional to the FN estimate: while at r a = 35 both the crystal 
and hybrid GWF tends to get delocalized, at r s = 60, the 
peaks of the hybrid GWF get more pronounced after the FN 
projection. 



uum (y — > 0) limit. Then, we determine precisely the 
critical value r* f» 31.5 ± 0.5 at which the transition be- 
tween the liquid and the hybrid phase takes place. After 
discussing the status of the mixed estimators, we end the 
section with showing that in addition to a upper bound 
to the energy, a lower bound can also be extracted from 
the FN-QMC datas. 

Comparaison with the litterature. We start with a dis- 
cussion of an (almost) square system of 56 particles in 
210 x 208 sites. This system is close to systems pre- 
viously studied in the litterature, and can hence serve 
as a reference. The absolute values of the FN ener- 
gies are given in Table. IVII together with the results of 
Ref£&. The same datas (without Ref£) are represented 
in Fig. 1111 where we plotted the zero-point motion en- 

1 /2 

ergy 2{Efn — Cir s /2]/r s (which converges toward C\ii 
at r s ^> 1) as a function of r s . We find a good agreement 
with the more recent datas of Refi& while the original 
datas of Ref£ are about ~ 0.015 higher (Ref£ and Ref£ 
used the same numerical code in their calculations). The 
observed small differences can be due to the following 
systematic errors: (i) too small < /3 sa t- Here we aver- 
aged the results over 50 < (3 < 100 and some more energy 
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r s = 20 


r s = 30 


r s = 35 


r 3 = 40 


r s = 50 


r s = 60 


r s = 70 


r 3 = 75 


Liquid 


56 


210x208 




-18.499 


-28.744 


-33.906 


-39.097 


-49.532 


-60.017 






Crystal 


56 


210x208 






-28.728 


-33.904 


-39.104 


-49.558 


-60.074 


-70.634 




Hybrid 


56 


210x208 


U* = 0.005 




-28.734 


-33.911 


-39.113 


-49.563 


-60.080 


-70.640 




Liquid 


57 


continuum 


Ref. j3 




-28.722 




-39.073 








-75.788 


Crystal 


56 


continuum 


Ref.T 




-28.700 




-39.090 


-49.526 






-75.904 


Liquid 


57 


continuum 


Ref.T 




-28.734 




-39.092 








-75.825 


Crystal 


56 


continuum 


Ref." 




-28.730 




-39.102 


-49.558 






-75.918 



TABLE I: Selected values of energy for the three GWF as well as values from the litterature for comparaison. The energies 
have been averaged over 50/t < < 100/t and calculated with an average of 25000 walkers. The statistical error is of the order 
of one in the last digit. The parameter A has been optimized at the variational level and is very well fitted by A » 0.59 r° a 57 
for the liquid and a little lower for the crystal while the liquid values have been used for the hybrid GWF. U* = 0.005 which 
corresponds to the splitting of the lowest band. The parameter do has also been optimized at the variational level and to very 
good approximation do/a ta (10/r s ) 1//4 . 



( 0.003) can probably be gained by using higher values of 
(3. Increasing (3 not only involves more computing time, 
but also using more walkers to keep the exponential in- 
crease of the variance under control, (ii) Small residual 
effect of the lattice (see Fig. ITU below ), (iii) shell effects 



(Refi&& used 57 particles (complete shell) for the liquid 
while we used the same N — 56 system for all three 
phases. This effect can account for a difference of energy 
of ~ 0.004 at r s = 40, see Fig. H21 below), and (iv) finite 
number of walkers. 
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FIG. 11: Comparaison with the results of the litterature for 
a system of 56 electrons in 210 x 208 sites. The plot shows the 
rescaled energy 2[£"pn — c\r s /2]/r^ 2 which converges toward 
the zero point motion C1/2 correction to the energy at high 
r s . The symbols correspond to our calculation with the liquid 
GWF (full circles), solid GWF (full squares), hybrid GWF 
(stars) and to the calculation of Ref. for the liquid (empty 
circles, 57 particles) and solid (empty squares). 



Thermodynamic (N — > oo ) limit. In Fig. 1121 we ex- 
amine the finite iV effect for the liquid GWF. Those 
can be already understood by looking at the finite N 
corrections to the energy at r s — (See Fig. ED Up- 
per Left): one finds that the (exact) kinetic energy 
E c (r s = 0, N) = 1 + AE C (N) strongly oscillates with N 
for small values of N before saturating toward 1 . These 
oscillations corresponds to the filling of the various shells 



of equal energies. Hence, following Ref£, we fit the N 
dependance of the liquid energy with the form, 

E(N) =E ao + E c (r s )AE c (N) - A(r s )/N (17) 

However, contrary to Ref. 0, we do not let the coeffi- 
cient E c (r s ) to be a fitting parameter but rather fix it 
to the calculated kinetic energy (for a given value of N, 
the finite N corrections to E c (r s ) are irrelevant at the 
precision of the calculations). Hence, we plot in the inset 
of Fig.^lthe liquid energy to which the E c (r s )AE c (N) 
correction has been substracted, and find a very good 
agreement with a 1 /N residual error. The quality of the 
overall (two parameters) fit can be found in the main 
Fig. E3 for up to 162 particles. In the Lowest Left panel 
of Fig. we have plotted the difference between the FN 
and the variational energy as a function of 1/N (circles 
for the liquid). We find that this difference presents some 
finite N effect about as strong as those of the FN energy 
alone (here £ FN « -39.088 + 2.5AE C (N) - 1.75/N while 
E VAR - E FN w 0.143 + 0.8A£ C (A0 + 1.3/iV.). The dif- 
ference between the finite iV behaviour of the variational 
and FN result is not very surprising since for those in- 
termediates values of r s (here r s = 40), the difference 
-Evar — -^fn is about ten time bigger than those effects. 
The same conclusion also holds for the hybrid (diamond) 
and crystal (squares) GWF. 

The finite N behaviour of the crystal and hybrid results 
are much simpler. They do not present the same oscilla- 
tory behaviour but rather converge smoothly toward the 
N ^> 1 value (see Upper Right panel of Fig. El) A good 
fit is obtained by plotting the FN energies as a function of 
1/N 3 / 2 (which is the expected finite effect for the crystal 
at large r s ). We find that the difference between crys- 
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FIG. 12: Finite N effect of on the Fixed Node energy with 
the liquid GWF at r s — 40. Energy averaged over 25 /t < f3 < 
bO/t as a function of the number N of particles with A = 4.5 
for a system of N = 2P 2 particles in 16P x 28P sites (P = 
2... 9 diamonds). The fit (solid line) corresponds to E = 
-39.088 + Ef{r s = 40)AP c (iV) - 1.7b/ N where Ef{r s = 
40) = 2.5 is the (fixed-node) kinetic energy for 50 particles 
and AE C (N) is the exact finite size correction to the kinetic 
energy computed at r s — 0. Inset: the oscillatory correction 
E^°(r s — 40)AE C (N) has been substracted from the energy 
which is plotted as a function of 1/N. The solid line is a linear 
fit to the data E = -39.088 - 1.75/7V. 



tal and hybrid slightly increases as N gets bigger (The 
respective fits are E FN = -39.108 - 10.1/7V 3 / 2 for the 
crystal GWF and E FN = -39.116 - 7.8/iV 3 / 2 for the 
hybrid.) We repeated the same procedure at r s = 50 
(up to 98 particles only) and arrive at a similar result 
(E FN = -49.56 - 13.8/iV 3 / 2 for the crystal GWF and 
E FN = -49.567- 12.2/iV 3 / 2 for the hybrid.) 

Continuum (v — > 0) limit. The presence of the un- 
derlying lattice can induce a correction to the contin- 
uum result. This effect is a small correction however and 
should not be mixed up with the much larger effects that 
can take plac o 32 i 33 at much higher values of r s . In addi- 
tion, contrary to the finite N effect which one would like 
to avoid (since we make use of Ewald resummation, wc 
do not look at mesoscopic samples in this study), there 
is an underlying lattice in real (semi-conductor based) 
samples. In Fig ll4l lower panel, we plot the energy as a 
function of the inverse of the surface of the sample for 
a fixed number of particle. We find indeed a l/(L x L y ) 
correction to the energy. However (Fig ll4l upper panel), 
this correction is almost GWF independant so that the 
relative stability of the phases is unaffected by the lat- 
tice for the v = 1/224 and v = 1/780 samples. For the 
latter, the correction is almost negligeable and we have a 
quantitative good agreement with the continuum model 
for the absolute values of the energy (see above, Fig. lllfl . 

Critical value of r s . We now determine precisely the 
critical value r* at which the liquid-hybrid transition 




FIG. 13: Finite N effect in a system of N = 2P 2 particles in 
16P x 28P sites. Upper Left panel: exact finite size correction 
AE C (N) to the liquid kinetic energy at r s — 0. This curve 
is used to remove the oscillatory dependance of the liquid 
energy due to the shell structure. Lower Left panel: difference 
between the variational anf fixed-node energy Pvar — Pfn as 
a function of 1/N for the liquid (circle, A = 4.5), hybrid 
(diamonds, U* — 0.12, A — 4.5) and crystal (square, do /a = 
0.74, A — 4). The linear fits (solid lines) are respectively y = 
0.143 + 1.3/iV, y = 0.13 + 1.3/N and y = 0.075 + 1.3/N. A 
better fit is obtained for the liquid using (dashed line) y = 
0.143 + 1.3/N + 0.8AE C (N). Right panels: FN energy for 
the crystal (squares) and hybrid (diamonds, U* = 0.12) as a 
function of -/V (Upper Right panel) and 1/N 3 ^ 2 (Lower Right 
panel). The fits (solid lines) are respectively y = —39.108 — 
10.1/7V 3/2 and y = -39.116 - 7.8/N 3/2 
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FIG. 14: Effect of dilution on the Fixed Node energy for a 
system of 30 particles in 40P x 42P sites (P = 1, 2, 3 and 4) at 
r s — 35. Lower panel: energy of the liquid GWF (circles) and 
solid GWF (squares, do = 0.82) as a function of the inverse 
volume l/(L x Ly). The lines are a guide to the eye. Upper 
panel: difference of energy between the liquid and the solid 
GWF as a function of 1/(L X L V ). 
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FIG. 15: Determination of the critical value r*. Left panel: 
energy of the liquid (full symbols) and hybrid (empty sym- 
bols) as a function of r s for f3 — 0,1,2,3 and 5/t (top to 
bottom) with their corresponding linear fits. The system con- 
tains 128 electrons in 128 x 224 sites and — cir s /2 + c 1/ / 2 ''s' /2 /2 
has been substracted to the energy. Upper Right panel: re- 
sulting r* (crossing points of the left panel) as a function of 
0. Lower Right panel rl (at — > oo) as a function of N for 
a filling factor of v = 1/56 (empty circles), v = 1/224 (full 
squares) and v — 1/780 (empty diamonds). 



takes place (Fig. 1151) . In the left panel of Fig.^J we plot 
the energy of the hybrid and liquid phase as a function 
of r s for various (rather small) values of (3. For each of 
them we extract the crossing point r* (/?) which is plotted 
in the upper right panel of Fig. We find that r* (0) 
saturates for rather small values of j3 ~ 2/t indicating 
that the gain in energy obtained by increasing (3 further 
(/3sat ~ 25 /t for this system) is almost GWF indepen- 
dant. The obtained r* is plotted in the lower part of the 
right panel as a function of the number of particles for 
several filling factors. At large N, r* converges toward 
r* = 31.5 ±0.5. 

Variational, mixed and fixed-node 
given observable O, three different 
constructed. The variational Ovar 



estimates. For a 
estimates can be 
(¥<dO|*G> 



and the the fixed node (forward walking) estimates 

O fn = (*™|0|*™) = (^le-^^Oe-^^l^c) 
have been used through this paper. An intermedi- 



FN\ 



ate one, the mixed estimate Omx = (^dO^g 
(^c\Oe~ H ™ P\tyc) is easy to compute and hence very 
common in QMC calculations. For most applications, 
\^g) is a reasonable approximation of the ground state 
I*" 



FN\ 



so that the mixed results Omx can be extrapo- 
lated toward the FN result using the (first order correc- 
tion) formula Ofn ~ 20mx~ Ovar- The mixed estimate 
is much easier to obtain than the FN estimate, since the 
latter requires to apply e~ H 13 ^ 2 after that the measure- 
ment has been made, and the deaths and births involved 
in the reconfiguration process considerably decrease the 
statistics. However, in the range of r s studied in this 




FIG. 16: Comparaison of the different estimates of the ki- 
netic and interaction energy for a system of 30 particles in 
40 x 42 sites at r s =25. Left panels: FN estimate for the liq- 
uid (diamonds), crystal (squares) and hybrid (circles) GWF. 
Right panels: mixed estimate for the liquid phase. Upper 
panels: kinetic energy E c (0). Lower panels: interaction en- 
ergy Eint(/3) 



work, we find important differences between the differ- 
ent estimates so that it is necessary to use the FN esti- 
mates. To illustrate this point, we plot in Fig. 11610pm (8) 
(Left panels) and Omx(P) (Right panels) for the kinetic 
energy E c (Upper panels, E c ~ 1 at r s = 0) and interac- 
tion energy E- lnt (Lower panels). For the case considered 
in Fig. 1161 the mixed estimate is a very poor estimate of 
the true (FN) estimate. Indeed, we have E^ AR ~ 2.01, 
Ef x {8 = 8) « 1.97 and E™(0 = 8) « 2.01 so that 
the interpolation from the mixed estimate would be very 
bad. 

Lower bound to the total energy. The evolution of the 
kinetic and interaction energy with 8 can be put to fur- 
ther use: we find that E^{8) < E^ h (f3) < Ef*{0) while 
E-^(0) < E^ t b (f3) < E^ t (0) (consistent with the inter- 
pretation that the hybrid phase is intermediate between 
the liquid and the crystal). As both £^ q (/3) and E^(8) 
increases with /3, (and since by construction the liquid 
and crystal GWF favor respectively the kinetic and inter- 
action energy), we conjecture that they form true lower 
bounds of the kinetic and interaction energies of the true 
ground state of the system. Hence, in addition to an up- 
per bound of the total energy, we can also construct a 
lower bound, 



£i"t y < E < E FN 



hyb 



(18) 



which for the case of Fig. Hfl gives -23.73 < E < -23.67 
with a precision of 0.3%. This lower bound is of interest 
for us, as it allows to estimate that the gain in energy 
provided by the hybrid GWF (with respect to the crystal 
and liquid one) is a significant fraction of the distance to 
the lower bond. 
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VII. DISCUSSION 

We conclude this paper with a discussion of several 
topics linked with the present study. 

Spin of the electrons. We restricted the present study 
to fully spin polarized electrons. Indeed, as the inter- 
action energy gets bigger, the system can minimize its 
interaction energy by antisymmetrysing the orbital part 
of its wave-function, and it was found (for instance in 
Refi&I) that above r s > 20 (which is the range of our 
study) the polarized liquid is more stable than the unpo- 
larized one. A strong in-plane magnetic field would also 
polarize the system at a smaller value of r s . 

Multiple exchange. At high r s when the crystal is well 
established, multiple exchange (of two, three or more 
particles) can take placed , leading to an effective in- 
teraction between the electronic spins. Recently these 
exchange interactions were estimated using the (finite 
temperature) path integral QMC techniqu e) 10 ! 35 . The au- 
thors proposed that the melting of the crystal is actually 
due to a divergence of multi-spin exchange. In those cal- 
culations, some sort of FN approximation was used be- 
low r s ss 60 to stabilize the crystal. Summing up the 
exchange energies calculated iniS, we find (in our unit) 
E cx as 0.001 at r s — 50 (where the calculation is most 
robust). This energy is relatively small and cannot ac- 
count for the difference observed between the hybrid and 
crystal GWF results. The difference between these path 
integral calculations and those presented in this paper 
lies in the FN approximation (the one used iniS seems to 
stabilize a well localized crystal) and the use of a rather 
large temperature T w 0.08 (which is roughly three times 
bigger than the energy difference between the crystal and 
the liquid at r s = 50 and a 100 times bigger than the cal- 
culated exchange energies). 

Hexatic and super solid phase. Among the proposed in- 
termediate phases is the presence of a hexatic quantum 
phase in between the liquid and the crystaliS^. Such 
a phase would have orientational order, but no transla- 
tional order. The hybrid GWF is made of Bloch waves 
that have, by construction, a hexagonal Fermi surface 
(one full band corresponds to momentum in the first 
Brillouin zone of the triangular crystal). However, we 
have also explicitely broken the translational symmetry 
so that, at the variational level, the hybrid GWF does 
certainly not correspond to a hexatic phase. Upon ap- 
plying the FN projection operator, the translational sym- 
metry tend to be restored (see the left panels of Fig. 
for instance) so that the present calculations cannot rule 
out the possibility of a quantum hexatic phase. Another 
proposition is the presence of a supersolid phaseSSiSi* 2 ^, 
i.e. the quantum coexistence of a liquid (of delocalized 
defects for instance) with the crystal. Qualitatively, the 
hybrid phase is in agreement with the idea of the su- 
persolid phase as discussed in Ref. |3|| In the supersolid 
described in Ref. however, one expects the number 
of electrons in the crystal to be smaller than the total 
number of electrons. We tried to construct such a GWF 



without success. 

Bubbles and Stripes. Around the melting point 
of a classical liquid-solid phase transition, the system 
gains some energy by being in a coexistence of the 
two phases. In the present system such a macroscopic 
phase separation is not possible since it would lead to 
a macroscopic dipole. Spivak and Kivclson recently 
proposed that a series of microemulsion phases (bub- 
bles/stripes of solid/liquid in the liquid/solid) can how- 
ever take placeiA£*i& If they exist, these microemulsions 
take place between the two values of densities n s and n\ 
(in electrons per bohr radius squared) that come out of 
the Maxwell construction: d n E sol (n s ) = d n E hq (ni) = 
[E liq (ni) - E so \n s )]/(ni - n s ) (energies in Rydberg). 
Solving for n s and n\ using the FN energies, we find 
that the (maximum) range of r s where these microemul- 
sions could occur is \r s — r*\ < Ar s with Ar s w 0.04 
(r* is the transition point. To a good approximation, 
Ar s w r* s oijc\ where the change of slope a w 0.003, de- 
fined as E llq w _E so1 + a(r s — r*) around the transition 
can be extracted from Fig|I]or Fig^J) The argument 
was originally developped for the liquid/crystal transition 
but holds equally well for the liquid/hybrid transition 
described in this paper (a direct liquid/crystal transition 
would lead to Ar s w 0.06). The maximum gain of energy 
given by the Maxwell construction (which uses the ener- 
gies calculated for neutral phases, hence only a fraction 
of this energy can possibly be gained by microemulsions) 
is AE « 2a 2 /(cir*) « 0.2/j.Ry (» 10~ 4 in our units). 

Metal-Insulator transition. About ten years ago, a new 
interest in two dimensional systems arised from the work 
of Kravchenko et at2i who reported on an unexpected 
metal-insulator transition. It was followed by an impor- 
tant body of literature but the origin of this transition 
is still under debate 3 ^ 3 ^ -. The problem involves under- 
standing the role of disorder in the system and lies outside 
the scope of this paper. Indeed the effect of disorder is a 
complicated issue and depends not only on the strength 
of the disorder but also on the typical length scale on 
which the potential varies. For instance, the component 
of the disorder varying on the scale of the inter electron 
distance is likely to pin and stabilize a crystal (or glassy) 
phase. A disordered potential varying on larger length 
scale however will favor fluctuations of the density and 
hence phase separation. In actual samples, long length 
scale fluctuations have been observed using local com- 
pressibility measurements* 1 -. 

A natural question for us is weither puddles of hy- 
brid/crystal phase are actually present in the experi- 
ments. As the mobility of the two dimensional systems 
improved in the past years, so had the critical value of 
r s where the transition is observed. In hole GaAs de- 
vices it goes from 41 i 42 r s = 20 — 25 to 43 ' 44 r s ~ 55 — 60 
in the cleanest samples, i.e. well above r* = 31.5. At 
r s = 60, the difference of energy between the liquid and 
the crystal phase is roughly 0.03 mRy per particle (see 
lower panel of Fig.^l so that the melting temperature of 
the crystal should be slightly larger (the crystal should be 
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stabilized by its spin entropy through the Pomeranchuk 
effect 20 ) which translates into temperatures of a few tens 
of mK. Those temperatures are therefore within the pos- 
sibilities of a good dilution fridge through slightly below 
the temperatures that are usually studied. 

Conclusion: how liquid is the hybrid phase? To sum- 
marize the main message of this article, we find that in 
these intermediate regions of densities, where little of the 
physics is known, the fixed node quantum Monte Carlo 
technique should be viewed as a probe. The freedom in 
the choice of the GWF allows one to start with a GWF 
that captures some sort of physics (as seen in the varia- 
tional calculation) and the evolution of the various (not 
only the energy) physical quantities upon applying the 
fixed-node projection operator gives us insights on the 
physics of the true ground state. In short, one introduces 
new physics in a GWF and sees weither this new physics 
is stabilized or washed out by the FN-QMC algorithm. 

With this paradigm in mind, we followed the evolu- 
tion (from variational to FN) of the physical quantities 



available within our algorithm. For 31.5 < r s < 60, we 
found that all of them show the same tendency to der- 
ealization. The succes of the hybrid GWF with respect to 
the crystal one is associated to its (much more) delocal- 
ized nature that allow a better balance between kinetic 
and electrostatic energy. The success of the hybrid GWF 
with respect to the liquid one however is associated to the 
change in symmetry as the hybrid GWF is constructed 
out of delocalized waves. We conclude that while the na- 
ture of the ground state above r* = 31.5 might still not 
be fully elucidated, the fact that it is not a liquid, yet 
not a localized crystal neither is now put on very firm 
grounds. At large r s the system eventually gets local- 
ized. At the present we cannot say weither this occurs 
through a crossover or a second transition. 
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